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ABSTRACT 



We investigate the value of horizontal turbulent diffusivity 7] by numerical 
calculation of thermal convection. In this study, we introduce a new method 
whereby the turbulent diffusivity is estimated by monitoring the time devel- 
opment of the passive scalar, which is initially distributed in a given Gaussian 
function with a spatial scale do. Our conclusions are as follows: (1) Assuming the 
relation r] — Lci'rms/S where v^uis is the RMS velocity, the characteristic length 
Lf. is restricted by the shortest one among the pressure (density) scale height and 
the region depth. (2) The value of turbulent diffusivity becomes greater with the 
larger initial distribution scale do. (3) The approximation of turbulent diffusion 
holds better when the ratio of the initial distribution scale do to the characteristic 
length Lc is larger. 

Subject headings: Sun: interior — Sun: dynamo — Stars: interiors 
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Introduction 



Turbulent diffusivity has been an important concept for tlie mean-field modeling of 



Miesch 



20051 ) ■ 



the interior convection and dynamo of the Sun and stars (see the review by 
It is a substantial factor for the transport of the angular momentum and magnetic field. 
While the non-turbulent molecular diffusivities are much smaller, i.e., molecular viscosity is 
z/ ~ 1 cm^ s~^ and molecular magnetic diffusivity is ?7 ~ lO'' cm^ s~^ in the solar convection 
zone, the random advective motion of gases in turbulence is considered to behave as a 
strong diffusion. The specific value is unknown but previous s tudies suggest that the value is 
around 10^°-10^^ cm^ s~^ 



le.g 



. Dikpati & Charbonneau 


1999 


). This value affects predictions 


(Choudhuri et al. 


2007; 


Dikpati & Oilman 


2006; 


Yeates et al. 



20081 ). It a lso affects the symrnetry o : 



polar field 



( iRempel 



Hotta fc Yokovama 



speed AQ (iHotta fc Yokoyama 



2010a 



the global magnetic field and the strength of the 



The value determines the difference in rotation 



20111 ) and the propagation speed of the torsional oscillation 



20071 ). Thus, the estimation of this value is crucially important. 



Some studies have alrea dy estimateci the y alue of turbulent diffusivity on the solar 



surface through observations. iWang et al 



( 119891 ) investi gated the ey o 



Chae et al. 



ution of active regions 



(120081 ) also estimated 



and derived the optimum value of turbulent diffusivity. 
the value of turbulent diffusivity through high resolution observations. They concluded that 
the turbulent diffu sivity depends on t he res olved scale, i.e., the value becomes smaller with 



higher resolution. 



Abramenko et al. 



( I2OIII ) also found this type of dependency through 



observation of bright points. 



Kapyla et al. 



( 120091 ) estimate the value of turb ulent diffusion with nu merical simulations 



of thermal convection using the test field method (jSchrinner et aL 



adopted for investigations of many different types of turbulence (IBrandenburg 



20051). which has been 



2005 



20081). 



A test magnetic field is passively transported by the convection fiows with no back reaction. 
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With the utilization of horizontally averaged values as a mean field, the coefficients of the 
g-effect and the tur bulent diffusivity are measured based on the mean-field equations. 



Kapyla et al 



( 120091 ) report that the value of turbulent diffusivity is proportional to the 



square of vertical velocity and is approximately proportional to the wavelength of the 



test ffeld. 



Cameron et al. 



(I2OIII ) investigate the value of turbulent diffusivity with a 



realistic radiative MHD simulation and estimate the value o 
the decreasing rate of the total magnetic ffux 



Yousef et al 



turbulent diffusivity from 



(I2OO3I ) estimate the turbulent 



magnetic diffusivity and kinetic viscosity in a forced isotropic turbulence using a similar 



way, i .e. from the decay rate of the magnetic field and the velocity field. 



Riidiger et al. 



(I2OIII ) use the cross helicity to e stimate the turbulent magnetic diffusivity in a stratified 



medium with forced turbulence. 



Riidiger et al. 



( I2OI2I ) extend this method to numerical 



calculation of thermal convection and observation of the sun. 

In this study, we introduce a new method to estimate the value of turbulent diffusivity. 
We investigate the development of a passive scalar whose initial condition is the Gaussian 
function. The method is found to be well suited for a Gaussian function at each time 
point and its peak density and spatial extent give us necessary information on the scalar's 
kinematics. A detailed explanation of the method is given in Section 12. 2[ The specific 
aims of this study are: (1) estimation of turbulent diffusivity of thermal convection with 
different sizes of the simulation box; (2) investigation into the validity of approximation of 
turbulent diffusion in thermal convection; (3) investigation into the dependence of turbulent 
diffusivity and the validity of approximation on the initial distribution scale. 
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Model 



2.1. Equations 



The three-dimensional hydrodynamic equation of continuity, equation of motion, 
equation of energy, and equation of state are solved in Cartesian coordinates (x, z) 
where x and y denote the horizontal directions and z denotes the vert ical direction. The 
formulations are almost the same as those used by 
expressed as. 



HottaetaL 



( 12012| ). Equations are 



^ = -V-[(po + Pi)v], 

_(v . V)v - - ^r/ez + -V ■ n, 

Po Po Po 

ds 1 

^ = _(v . V)(so + si) + — • (i^poToVsi) 

at pqTq 
pi 



7 



Po 



(n-v)-v. 



Pi = Po 7 \- si 

Po 



(1) 
(2) 

(3) 

(4) 



where Po(-2), Pol^^), Tq{,z), and Sq{z) denote the time-independent, plane-parallel reference 
density, pressure, temperature, and entropy, respectively and denotes the unit vector 
along the z-direction. 7 is the ratio of specific heats, with the value for an ideal gas being 
7 = 5/3. pi, pi, and Si denote the fluctuations of density, pressure and entropy from 
reference atmosphere, respectively. Note that the entropy is normalized by specific heat 
capacity at constant volume c^. The quantity g is the gravitational acceleration, which is 
assumed to be constant. The quantity 11 denotes the viscous stress tensor. 



Po^^ 



dvi dv 



dxj dxi 3 

and V and K denote the viscosity and thermal diffusivity, respectively, v and K are 
assumed to be constant throughout the simulation domain. 



(5) 



We assume an adiabatically stratified polytrope for the reference atmosphere except 
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for entropy: 



Ho{z) 



Pr 



T 



(m + l)Hr 

z 

(m + l)Hr 

z 

(m + l)Hr 



m+l 



pog' 



(6) 

(7) 

(8) 
(9) 



where pr, Pr, T^., and Hj. denote the values of po, po, To, Hq (the pressure scale height) at the 
bottom boundary z = 0. The profile of the reference entropy so{z) is defined with a steady 
state solution of the thermal diffusion equation V • {K pqTq'V sq) = with constant K: 

7(5(z) 



dso 
dz 

5{z) = S, 



Ho{z) 

pr 

>o(^) 



(10) 

(11) 



where S is the non-dimensional superadiabaticity and 5r is the value of 5 at 2; = 0. In 
spite of a non-zero value of superadiabaticity, the adiabatic stratification is acceptable due 
to the small value of superadiabaticitiy. The strength of the diffusive coefficients u and 
K are expressed with the following non-dimensional parameters: the Reynolds number 
Re = VcH^/u, and the Prandtl number Pr = I'/K, where the velocity scale = 
In all cases of this study, the parameters are set as Re = 300, Pr = 1, 5j. = 1 x 10~^. We 
calculate three cases with different box sizes (see Tabled]). The horizontal size is the same 
in all calculations, i.e. = Ly = L = 26.16ifi. and the number of grids in x, y directions 
are set as N^. = Ny = 1152. We adopt three different vertical sizes of box, = 2.18ifr, 
1.635ifr, and 1.09ifr for cases 1, 2, and 3 respectively. The number of grids in these cases 
are set as = 96, 72, and 48, respectively. The Rayleigh number, which is defined as 



(12) 
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in these cases are estimated to be 1.3 x 10^, 3.4 x 10^, and 1.7 x 10^ respectively, where 
As denotes the difference of entropy between the top and the bottom boundaries. The 
calculation domain is —L/2 < x < L/2, —L/2 < y < L/2 and < z < L^. The bounda ry 



HottaetaL 



torn . 



conditions and the numerical method are the same as those used by 
The boundary condition for the x, y direction is periodic for all variables, and the stress free 
and impenetrative boundary conditions are adopted and the entropy is fixed, i.e. Si = at 
z = and L^. 



2.2. Method for Estimation of Turbulent Diffusivity 

In this study, we calculate the evolution of passive scalar to estimate the value 
of turbulent diffusivity. Along with the equations ([I])-(|4]), we simultaneously solve the 
advection equation of the passive scalar as 

^ = -v-(gv), (13) 

where Q is passive scalar density. Although in eq. (fT3|) . the diffusion term does not appear 
ex plicitly, we use tiny a rtificial viscosity on the passive scalar, a technique which is adopted 



in (IRempel et al 



20091). Its initial condition is set as 

2 



2 2 

Q(x, y,z,t = 0) = exp ( ) . (14) 



We adopt three initial conditions, i.e. do = 2.5H,., 5.0H,., and 7.5H,. for each of the 
different depth settings (cases 1-3); hence the total number of cases is nine. In the initial 
condition the passive scalar does not depend on z, since we focus on the turbulent diffusion 
in the horizontal direction. Since the transport of the passive scalar is assumed to be 
approximated by a diffusion process with constant diffusivity rj, then its density should 
obey the two-dimensional diffusion equation as 

dQ fd'd'\ 



dt \ dx"^ dy"^ 
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When the calculation domain is infinite, the analytical solution of eq. f lTSj) with the initial 
condition of eq. ([H]) is expressed as 

O-f^Vexpf-^V (16) 



where, d'^ = Arjt + dl. In this study we adopt periodic boundary conditions; thus the analytic 
solution is given by the periodic superposition of the above formula and can be expressed as 



«= E E 



exp 



d^ 



(17) 



i=— oo j=— oo 

When the width of the Gaussian function is narrower than box size {d < L), the analytical 
solution in the range, —L/2 < x < L/2 and — L/2 < y < L/2, can be approximated as 



i=-ii=-i 



do 

■ exp 



rf2 



We estimate the value of turbulent diffusivity by the following steps: 



1. The advection eq. f lT3|) is calculated with the obtained velocity of thermal convection. 

2. The obtained passive scalar in each step is vertically averaged as 



1 



Q = -r Q^^- (19) 

Jo 

Note that by using this method, we will obtain an averaged turbulent diffusivity along 
the 2;-direction. 

3. The result of averaging, i.e., eq. ( !T9|) . is fitted with eq. (fT8|) . Note that the fitting has 
only one parameter d{t), and this parameter has information on both the height and 
the width of the Gaussian function. 

4. According to the analytical relation, rf^ = Arjt + cig, we obtain the value of turbulent 
diffusivity from the slope of d'^{t). 
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3. Results & Discussion 

Figure [1] shows the results of our hydro dynamic calculation. The three panels in the 
left, middle and right columns show the contours of entropy in cases 1, 2 and 3, respectively. 
Due to the large Rayleigh number, the velocity with a large box size is high (see the third 
row in Table [1]). A detailed investigation of cell size distribution will be reported in our 
forthcoming paper (lida et al, in preparation). 

Figure [2] shows the contour of the passive scalar whose width of the Gaussian function 
at the initial condition is = 2.5ifr- We can see that the passive scalar is diffused with 
turbulent convection. The dependences of d'^ on t are provided in Figure [3l and are shown 
to be almost linear. This shows the validity of the diffusive description for the turbulent 
transport by the convective motion. The estimated turbulent diffusivity is shown in Figure 
H^. It is derived through linear fittings to the curves in Figure |3] in the range of < t < tmax 
where tmax is chosen so as to reduce the fitting error; it is given in Table [TJ 

The scaling behavior of the obtained diffusion is studied by changing the depth of 
the simulation box in cases 1, 2 and 3. In Figure UK, the blue, green, and red lines show 
the values of turbulent diffusivity with do = 7.5Hj., 5.0H^, and 2.5Hj., respectively. The 
value of turbulent diffusivity scales with the size of the box, which is discussed in the 
next paragraph. The value of turbulent diffusivity also scales with the initial width of 
the Gaussian function do. Using a wider Gaussian function makes the larger size of the 
convection cell work more efficiently and generates a larger value of turbulent diffusivity. 

In the mean field model, it is thought that the coefficient of turbulent diffusion can 
be expressed as r/ = Lcfrms/3, where Lc is the characteristic length scale of turbulence and 
frms is the root-mean-square (RMS) velocity. The value of turbulent diffusivity is obtained 
in this study, and we can estimate the value of based on the mean field model i.e., 
Lc = 377/frms- The estimated horizontal RMS velocities and characteristic length scale are 
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shown in Figure |Dd and c, respectively. We discuss the dependence of Lc on the box size 
with (io = S.Oifr and 7.5H^. With a smaller box, i.e. l.GSbH,. (case 2) and l.OQifr (case 3), 
the characteristic lengths are almost the same as the sizes of the boxes (the size of the box 
is indicated by the dashed line). It is natural that the largest cell size is determined by the 
size of the box and that the largest cell is most effective for advecting the passive scalar. 
Although we expected that the characteristic length of case 1 would also be the same as 
Lz, the obtained characteristic length scale was smaller than even with do = S.Oifr and 
l.bHj.. A possible reason for this result is that the characteristic length Lc is restricted by 
the convection cell, which is also limited vertically by the pressure scale height {H^) or the 
density scale height {'^H^ ~ LQIH^). Although should be evaluated in the horizontal 
scale, the mixture of the passive scalar may occur at approximately the same distance with 
the vertical scale. It should also be noted that when the narrowest Gaussian function, i.e. 
(io = 2.5Hj. (red line), is used, the characteristic lengths are restricted by the width of the 
Gaussian function in cases 1 and 2. 

Next, we discuss the validity of the approximation of turbulent diffusivity quantitatively. 
We calculate the estimated error of the linear fitting of (P as: 



where N is the number of data points along the time, d?{tn) is the n-th estimated result 
and d?{tn) is the n-th result of the fitted line. In Figure Htl, we found a dependence of a 
on (io, i.e. a is larger with narrower d^. Although the qualitative relation is not clear, it 
indicates that the estimated error a tends to become smaller with a larger ratio d^/Lc of 
the initial width of Gaussian function to the characteristic length (Lc) 




(20) 
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4. Summary 

We investigated the value of horizontal turbulent diffusivity r/ by a numerical calculation 
of thermal convection. In this study, we have introduced a new method, whereby the 
turbulent diffusivity is estimated by monitoring the time development of the passive scalar, 
which is initially distributed in a given Gaussian function with a spatial scale do. Our 
conclusions are as follows: (1) Assuming the relation t] = Lcfrms/3 where frms is the RMS 
velocity, the characteristic length is restricted by the shortest one among the pressure 
(density) scale height and the region depth. (2) The value of turbulent diffusivity becomes 
larger with a larger initial distribution scale do. (3) The approximation of turbulent 
diffusion holds better when the ratio of the initial distribution scale do to the characteristic 
length Lc is larger. 



Conclusion (2) is c onsistent with the results of obse rvational study (jChae et al.ll2008l : 



Abramenko et al. 



20091). In this study. 



20111 ) and a previous numerical study (iKapyla et al. 
we do not estimate the correlation length directly from the thermal convection. This will be 
achieved in our future work with an auto-detection technique and our characteristic length 
{Lc) will be compared with directly estimated correlation length. We now assume that 
our characteristic length is an average of correlation length at each height (lida et al. in 
prepareation) . The turbulent diffusion in the horizontal directions is estimated in this work, 
but such estimations are also important in the vertical directions for addressing the solar 
dynamo problem from the viewpoint of the transport of magnetic flux from the surface to 
the bottom of the convection zone. Such a study will be conducted in the future. 

Although turbulent diffusivity averaged in the whole box is estimated in this study, the 
dependence of this estimation on the height is important. There are, however, two reasons 
why it is difficult to estimate this dependence with our method. First, in our calculations 
the integrated passive scalar density is not conserved at each height. Second, we found that 



it is difficult to estimate the diffusivity separately for each horizontal plane only by solving 
eq. f[T^ two-dimensionally in the x — y plane the because the results show that the passive 
scalar density is strongly concentrated in the boundaries of the convection cells. Such a 
spatially intermittent structure is inappropriate for obtaining a statistical property like 
the turbulent diffusivity. These difficulties will necessitate some substantial improvements 
in our method. We are also interested in the effect of feedback from th e magnetic field to 



the convection beca use of its influence on the turbulent diffusivity (e.g. 



Riidiger et al. 



Yousef et al. 



2003 



20111 ). 
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Table 1: Parameters of study's calculation. 



Case 1 



L^xLyX (H^) 26.16^ x 2.18 26.16^ x 1.635 26.16^ x 1.09 

N^xNyX 1152^ X 96 1152^ x 72 1152^ x 48 

Po{L,)/poiO) 22 4.9 2.4 

^^rms (^^c) 0.287 0.206 0.148 

W {Hjv^) 75 112.5 150 
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Fig. 1. — Contour of entropy (si). The three panels in the left, middle, and right columns 
correspond to the results of cases 1, 2, and 3, respectively. The rows in the top, middle and 
bottom columns show the plot at z = Lz, Lz/2, and 0, respectively. 
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Fig. 2. — Contour of passive scalar at t = 75Hj./vc- The panels in the left, middle, and right 
columns show the results in cases 1, 2, and 3, respectively. The panels in the top row show 
the contour of passive scalar Q at z = 0. The bottom row shows the plot of passive scalar 
density averaged over z, i.e.,Q defined by eq. ffT^ . 



(a) L,=2.18H, 




(b) L,=1.635H, 
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c) L,=1.09H, 
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Fig. 3. — Distribution range cP of the passive scalar as functions of time. Panels a, b, and 
c show the dependence of rf^ with the sizes of boxes = 2.18ifr, 1.635ifr, and 1.09H,., 
respectively. The blue, green, and red lines show the results with do = 7.5Hr, S.Oifr, and 
2.5Hj., respectively. 
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(a) Turbulent diffusivity 






(c 




2.4 : 




2.2 : 




2.0 : 




1.8 : 


o 

_l 


1.6 ■ 




1.4 ■ 




1-2 ■ 




i.n 



1.0 1.2 1.4 1.6 1.8 2.0 2.2 

Lz (H.) 

(c) Characteristic lengtli 




1.0 1.2 1.4 1.6 1.8 2.0 2.2 

Lz (H.) 



(b) Horizontal RMS velocity 

0.35 




1.0 1.2 1.4 1.6 1.8 2.0 2.2 

Lz (H,) 
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Fig. 4. — (a) Dependence of turbulent diffusivity on the size of box. The value of turbulent 
diffusivity is normalized by rjo = v^H^ (b) Dependence of horizontal RMS velocity on the 
size of box. (c) Dependence of characteristic length on size of box. The dashed line shows 
the size of box. (d) Dependence of estimated error on the width of the Gaussian function. 
In panels a and b, the blue, green, and red lines show the results with rfo = 2.5H^, 5.0H,., 
and 7.5Ht-, respectively. In panel d, the solid, dotted and dashed lines show the results with 
= 2.18ifr, 1.635ifr, and l.OQi^i., respectively. 



